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The localized Hartree-Fock potential has proven to be a computationally efficient alternative to the optimized 
effective potential, preserving the numerical accuracy of the latter and respecting the exact properties of being 
self-interaction free and having the correct —1/r asymptotics. In this paper we extend the localized Hartree- 
Fock potential to fractional particle numbers and observe that it yields derivative discontinuities in the energy 
as required by the exact theory. The discontinuities are numerically close to those of the computationally 
more demanding Hartree-Fock method. Our potential enjoys a “direct-energy” property, whereby the energy 
of the system is given by the sum of the single-particle eigenvalues multiplied by the corresponding occupation 
numbers. The discontinuities c-j- and cj, of the spin-components of the potential at integer particle numbers iVj- 
and IV j, satisfy the condition cqlVj- + c4.IV4 = 0. Thus, joining the family of effective potentials which support 
a derivative discontinuity, but being considerably easier to implement, the localized Hartree-Fock potential 
becomes a powerful tool in the broad area of applications in which the fundamental gap is an issue. 
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I. INTRODUCTION 

The total energy of a quantum-mechanical system, re¬ 
garded as a function of the continuously varying num¬ 
ber of electrons, is a series of straight line segments 
with abruptly changing slopes at the integral values of 
the particles number 1 . Although this function is con¬ 
tinuous, its derivative has jumps at integral values of 
the number of particles - a property that is usually 
referred to as the derivative discontinuity. In density- 
functional theory (DFT) 2,3 the derivative discontinuity 
in the energy is complemented with a discontinuity in 
the exchange-correlation (xc) potential v xc (r): this po¬ 
tential experiences a jump when the number of particles 
passes through an integral value 4 . 

The role of the derivative discontinuity in DFT is crit¬ 
ical for the correct interpretation of the fundamental 
gap 4,5 and in quantum transport 6,7 , to name only two 
areas. Simple approximate DFT schemes relying on lo¬ 
cal or semi-local xc functionals do not, however, satisfy 
the derivative discontinuity requirement 1 . On the con¬ 
trary, the optimized effective potential (OEP) 8,9 stands 
out as the exact exchange potential which generally per¬ 
forms very well both in DFT and in time-dependent DFT 
(TDDFT) 10-12 and, in particular, supports the derivative 
discontinuity 13,14 . In principle, OEP provides a natu¬ 
ral starting point for further systematic inclusion of the 
correlations within DFT 15 . The notorious drawback of 
OEP is, however, the extremely high computational cost 
of its implementations, which is due to the necessity to 
know both empty and occupied orbitals and solve the so- 
called OEP integral equation. Approximations to the ex¬ 
act OEP equations 16,17 reduce the numerical effort, but 
can hardly serve as a solid basis for further development 
of the theory, either in the direction of including correla¬ 
tions, or for dealing with time-dependent phenomena. 

An attractive alternative to OEP, requiring the knowl¬ 
edge of occupied orbitals only, was proposed by Della 


Sala and Goer ling 18 , and it is known to yield accu¬ 
rate results for both closed and open shell atoms and 
molecules 18-20 . This potential is generally known in the 
literature as localized Hartree-Fock (LHF) potential. Re¬ 
cently, the LHF potential has been extended to TDDFT 
and has been successfully used to compute the dynamic 
response of the interacting electron gas 21 . Similar to the 
OEP, the LHF potential satisfies important requirements 
of the exact DFT: it is self-interaction free, it has the cor¬ 
rect —1/r asymptotic behavior, and, as recently shown, 
it can be rigorously derived from a minimum variational 
principle within the optimized-propagation scheme 21 . In 
this paper, one more fundamental property is added to 
this list: we report the success of the LHF potential in 
dealing with fractional particle numbers in atomic sys¬ 
tems, producing derivative discontinuities comparable to 
the Hartree-Fock (HF) method. 

We also point out some subtle but conceptually im¬ 
portant differences between the LHF potential and a 
standard DFT potential, by which we mean a poten¬ 
tial that can be expressed as a functional derivative of 
an exchange-correlation energy functional. For spinless 
fermions, or, equivalently, for fully spin-polarized sys¬ 
tems, the LHF potential is uniquely determined by the 
variational principle and thus suffers no discontinuity, 
even though the energy has a derivative discontinuity. 
This is possible because the LHF potential is not the 
functional derivative of the energy functional. In fact, 
we show that the LHF is a “direct-energy potential”, in 
the sense of yielding the energy as the sum of the single¬ 
particle eigenvalues multiplied by the corresponding oc¬ 
cupation numbers 22 . In the general case, when both spin 
components are present, we find that the LHF potential 
experiences a discontinuity whenever the particle num¬ 
bers in both spin components are integers. This is similar 
to DFT, but again there is a difference: the up-spin and 
down-spin components of the DFT potential at integer 
particle numbers are defined up to two arbitrary inde- 
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pendent constants, whereas for the LHF potential the 
two constants are linked by a constraint, which reduces 
the degree of arbitrariness. 

This paper is organized as follows. In Sec. II, we 
summarize previous results of the optimized propagation 
method for integral number of particles and we extend 
this method to a statistical mixture with fractional oc¬ 
cupation numbers, emphasizing the general properties of 
the LHF potentials. We also show that the total en¬ 
ergy is simply given by the sum of the LHF eigenvalues 
multiplied by the corresponding occupation numbers. In 
Section III we present the results of analytic calculations 
performed with the LHF potential with fractional par¬ 
ticles number in certain simple cases. In Sec. IV, we 
present and discuss results of numerical solutions per¬ 
formed with the LHF potential for atomic systems with 
fractional particles number. Section V contains conclu¬ 
sions, and the proofs of some technical facts and detailed 
derivations are collected in the Appendices. We use the 
atomic units (e 2 = H = m e = 1) throughout. 

II. FORMALISM 

First we summarize the essentials of the LHF method 
for integral particle number TV. Our system is described 
by the many-body Hamiltonian 


Hn = T n + V ext ,N + Un, 


where 



*=l 


is the kinetic energy operator, 


N 

V e xt,N = ^ ' ^ext{^i) 
i=1 


is the external potential energy operator, and 


U N = 


N 

£■ 

i<j 


( 1 ) 

( 2 ) 


( 3 ) 


( 4 ) 


is the Coulomb interaction energy operator. 

The localized Hartree-Fock (LHF) potential, initially 
introduced by Della Sala and Gorling 18 , has recently 
been re-thought in the more general context of time 
propagation, as the single-particle potential whose time- 
dependent Slater determinantal solution comes closest 
to fulfilling the many-body time-dependent Schrodinger 
equation 21 . Equation (6) of Ref. 21, which we take as 
the starting point of the present work, is a self-consistent 
equation for the LHF potential, represented in the form 


v ef f(x) = v ext (x) + v(x), (5) 


with x = (r, cr) standing for both space and spin coor¬ 
dinates, and where the potential v(x) is found from the 
self-consistent equation 

~t;(®n\pn(x)(V n — U N )\$ N ) = 0 , (6) 


where 


N 

V Ar = ^fi(x i ), (7) 

i =1 

N 

pN (x) = ^ S ( r i “ r ) ,<r (8) 

2=1 

is the spin-resolved density operator, and \$n) is the 
ground state of the effective noninteracting Hamiltonian 

H-ef f,N = T N ~b Vext,N + Vn- (9) 

Notice that there are two equations in Eq. (6), one for 
each spin orientation, which determine the two compo¬ 
nents of the LHF potential, fi(r,f) and fi(r,),). The 
physical content of these equations, in addition to the 
points discussed in Refs. 18 and 21, is that the expec¬ 
tation value of Vn coincides with the expectation value 
of the Coulomb interaction when evaluated on the sub¬ 
set of configurations that have one particle of spin a at 
position r the probability of each configuration being 
determined by the wave function of the noninteracting 
ground state |$jv). Since \$n) is a single determinan¬ 
tal state, the expectation value of the two-body operator 
Pn(x)Vn and the three-body operator pn(x)Un can be 
straightforwardly evaluated by means of Wick’s theorem, 
and expressed in terms of the average spin resolved den¬ 
sity njv(x) and the density-matrix pvr(x, x') of the nonin¬ 
teracting ground state, leading to explicit self-consistent 
equations for v(x), as shown in Ref. 21. 

In order to extend the formulation to fractional parti¬ 
cle numbers TV + a, with 0 < a < 1, we introduce the 
weighted average of Eq. (6), namely 

— Tf — {*&n\pn (x) (Vn - &n)\&n)+ 

£ (10) 

-^^{4w+i|/5at+i(x)(V)v + i - Hjv+OI'lW+i) = 0 , 

where |3 >jv) an d |3 >jv+i) are the ground states of the effec¬ 
tive noninteracting Hamiltonian H e f f with N and N + 1 
particles, respectively. Notice that all orbitals are calcu¬ 
lated with the same potential v(x), and |$jv) and |$jv+i) 
differ only in the TV+1-th orbital, which is empty in |$jv) 
and occupied in |$jv+i)- This orbital has a definite spin 
orientation: therefore our formulation can accommodate 
fractional occupation of one spin component or the other, 
but not of both simultaneously. (Note: the above formula 
is valid for TV > 1. For TV = 0 only the second term is 
present) 
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Equation (10) is very similar to what one would obtain 
by replacing the average in I^at) in Eq. (6) by the average 
in the fractional ensemble Dn+<x = (1 — a)|<frjy)(<&jv| + 
a|$jv+i) (^jv+i |- There is a subtle difference, however, 
arising from the choice of different normalization factors 
N and TV + 1 in the two terms of Eq. (10). Besides 
following naturally from the time-dependent formulation 
of Ref. 21, this choice guarantees that, upon integrating 
Eq. (10) over x (i.e., integrating over r and summing up 
over a) we get 


(1 — a)($jv|Vjv — £/jv|3bv) + 
aK^iv+il^v+i — E/jv+iI^jv+i) = 0, 



a 

FIG. 1. The renormalized occupation fraction of HOMO /3 
[Eq. (13)] used in the calculation of the effective potential 
versus the physical fraction a. 


that is to say, the average of V in the fractional ensemble 
Dn +0 coincides with the average of U in the same en¬ 
semble. This means that V is a “direct-energy” potential 
in the sense of the recent paper by Levy and Zaharias 22 , 
i.e., in the sense that the expectation value of the many- 
body Hamiltonian H in the ensemble D^ +a will be given 
by the sum of the single-particle eigenvalues of H e f f mul¬ 
tiplied by the corresponding occupation numbers 


N 

\ " 


En+o. — 7 y G + OiCN+1- 

i= 1 


( 12 ) 


The slight complication caused by the normalization fac¬ 


J 


[^ +/3 ( x ) + Gn+p] n w+/ 3 (x) = / 


tors in Eq. (10) can be easily circumvented by noting that 
the left hand side of that equation is proportional to the 
expectation value of p(n)(V — U) in the fractional ensem¬ 
ble b N+ p = (1 - /?)|$jv)($jv| + /3|$AT+i)($Ar+i| where 


P = 


aN 


(13) 


(1 — a)(l + N) + aN ' 

Then, our generalized LHF equation takes the final form 

(p(x)(V-U)) N+ p = 0, (14) 

with /3 given by Eq. (13). In Fig. 1 we plot the renor¬ 
malized occupation fraction, /?, versus the physical one, 
a, for several values of N. 

The evaluation of the expectation values of the two- 
and three-body operators that appear in Eq. (14) is 
greatly simplified by the use of a generalized Wick’s the¬ 
orem, which we prove in the Appendix B. The theorem 
states that the average of a product of field operators in 
a fractional ensemble such as Dn+p can be calculated as 
a product of averages of pairs of field operators in the 
same ensemble. The net result, derived in details in the 
Appendix A, consists of two equations, (15) and (16), 
which must be satisfied simultaneously 


„ N +t 3 


(*')- 


r - r' 


|p Ar +/ 3 (x,x')| 2 dx' 


+ 


p^(x, x'V^+V, x")p jv+ ^(> 


(15) 


-dx'dx", 


G N+a = 0. (16) 

In Eqs. (15)-(16), 

N 

p N+ P(x,x.') = ^<7i(x)0-( x/ ) + ^JV+l(x)^ + i(x') 

2=1 


( 17 ) 
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is the density-matrix, where 0,(x) are spin-orbitals, 


n JV+/3( x ) = p N +Pfa X ) 


(18) 


is the spin-resolved density, 


,f+0( x )=i»+0(x)-»W(r) 


(19) 


is the exchange potential, where n^ +/3 (r) = J dx.' n ^ is the Hartree potential, and 

y (x)n JV+7 (x)dx+i J ^ +7 (r)n^ +7 (r)dr+i J l^ +7 ( x ’ x ')| 2 


Gat+7 = / ^ +7 ^ 


dxdx' 


r - r' 


( 20 ) 


is a 7-dependent constant. It is also shown in the Ap¬ 
pendix A that in all cases, except for the fully spin- 
polarized one, 


and n Ar+/3 (x) all the orbitals except for the HOMO of 
each spin orientation. Substituting into Eq. (15), we find 
straightforwardly 


G n+p = 0, (21) 

which both simplifies Eq. (15) and can be conveniently 
used instead of Eq. (16). 

It is important to notice, at this point, that our equa¬ 
tions determine v x uniquely only when the particle num¬ 
ber is fractional. Otherwise, when both Nf and ZVj. are 
integers, the transformation 

^(r,t) —> v*(r,t) + Cf, 

Uz(r,4.) -> Ua.(r,4.) + ci, 

where c-j- and cj. are constants, leaves both equations sat¬ 
isfied if the condition 


c-f-ZV-j- + cj.1V.j, — 0, (23) 

is fulfilled. (This can be verified explicitly using the idem- 
potency of the density matrix p(x, x') for integer parti¬ 
cle number). Therefore, at integral and only at integral 
number of particles and if both N- j- and IVj. are non-zero, 
the exchange potentials v x (r, a) are not defined uniquely, 
and we will see below that they experience a jump when 
the particle number passes through an integer value. 

In the special case of a fully spin-polarized system (say, 
= N is a non-zero integer and IVj. = 0), w x (r,j') 
is uniquely defined, while ^(r, 4.) is not defined at all 
and is irrelevant. Letting the number of particles tend 
to N from below and above, while maintaining the full 
spin-polarization, we see that Eqs. (15) and (16) tend to 
the same equations with N particles, the latter having a 
unique solution v x (r, f). Therefore, in this case the ex¬ 
change potential does not have a jump at integer particle 
number. We will see in Sec. IIIC that there still is a 
jump of the total energy derivative in this case, leading 
to the important conclusion that the discontinuity of the 
energy derivative and that of the exchange potential are 
not directly related properties (cf. Ref. 22). 

To determine the asymptotic behavior of the poten¬ 
tial far outside the system, we neglect in p N+,3 (x, xi) 


^ +/3 (r,cri)-^ +c ai , 

v x +/3 (r,a 2 ) -^ +c CT2 , 


(24) 


where 07 is the spin of the fractionally occupied HOMO, 
u 2 is the opposite direction of spin, and c CTl 2 are con¬ 
stants. 

While within our approach the energy of a system is a 
sum of the orbital energies multiplied by the occupation 
numbers, the standard DFT expression for the energy 
holds as well. That expression reads 


N 

En+cx = 'y ] Ci + aCN+1 

i=l 


v^ +a (x)n N+a (x) dx 


1 

2 


v H +a ( r )n N+a (r)dr + E?+ a , 


(25) 


with 


E N+a = _ 


1 /• |p A '+ a (x,x' 


,'M2 


r - r' 


-dxdx', (26) 


while the sum of the last three terms in Eq. (25) is zero 

E X **~J'v^° (x)n N+a (x)dx -I/) (r)n N+a {r)dr = 0, 

(27) 


by virtue of Eq. (16). 

Recently, in the framework of the general DFT, it has 
been shown that, by shifting the exchange-correlation po¬ 
tential by an appropriate constant, it is always possible 
to make the total energy equal to a sum of the orbital 
energies 22 . In our case, we have the latter property au¬ 
tomatically built in the formalism. Furthermore, follow¬ 
ing the line of arguing in Ref. 22, we note that the or¬ 
bitals, and hence the density n(x) and the density-matrix 
p(x, xi), must be continuous across the integer number of 
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particles. Therefore, by Eqs. (26) and (27), we conclude 
that the quantity 

J v x (r)n(r)dr= J v x (r, f)n(r, f)dr + J v x (r, |)n(r, |)dr 

(28) 

is continuous with respect to the particles number. Then 

/A« x (r, t )n(r,t)dr + J Av x (r, l)n(r, i)dr = 0, (29) 

where Av x are the jumps in the potentials when the parti¬ 
cle number passes an integer value. Since these jumps are 
independent of r where the corresponding spin-densities 
are non-vanishing, we can rewrite Eq. (29) as 


A v x ^N^ -(- A v x iN± = 0. 


(30) 


Eqs. (15) and (21) can be solved to 24 


zt( r ) = - J 


i) 


|ri - r 

<Ai( r i) 


dr i + c^, 


v xi (r) = -p [ - dri, 

J l r i — r | 


ct = -ej 


<t> t( r )^( r x) 

lri — r| 


drdr\ 


a 


2~a 


Since 

vh( r) = 

we also have 


+ P4>\{ r i) 


l r i — r 


dr\, 


(34) 

(35) 

(36) 

(37) 

(38) 


From Eq. (30), in the fully spin-polarized (“spinless”) 
case (say, = 0) we have the continuity of the exchange 
potential (At; x <|- = 0) with respect to the variation of 
the particles number. Otherwise, if both Nf and TVj, are 
non-zero, there can be discontinuities in the exchange 
potentials, while Eq. (30) must be satisfied. Below we will 
see examples of the realization of the both possibilities 
(Figs. 2 and 5, respectively). 

Recently an approach to the derivative-discontinuity 
problem was proposed which takes use of the ensemble 
generalization of an arbitrary DFT xc functional 23 . We 
note that our method is within the same lines, being the 
ensemble generalization of the LHF theory. 


«t( r ) = P J |^_ 1 y| rfr i + c t> ( 39 ) 

M r ) = f ( 4 °) 

By Eq. (12), we can write for the total energy 

E = e t + CMq. (41) 

Let a (and, consequently, /3) be small. By the perturba¬ 
tion theory to the first order in a we can write 

E = e | + ae°, (42) 


III. ANALYTICALLY SOLVABLE CASES 

A. Number of particles between 0 and 1 

To consider the case of a particles, 0 < a < 1, we must 
use Eq. (10) with only the second term. Therefore, we 
have 

v(r) = 0, (31) 

leading to 

«x(r) = -vii(r), (32) 

Veff(r) = v ext (r), (33) 

all of which are exact results. 


where e2 and ej are the eigenvalues of the Hamiltonians 

- ^A + v ext (r), (43) 

- ^A + v ext (r) + J ^°^| dri, (44) 

respectively, and </>o( r ) is the ground-state eigenfunction 
of a particle in the potential v ex t(r). Equation (42) is 
obtained with the use of Eqs. (39) and (40) noting that 
by Eq. (36) the first-order term in a vanishes for e-j-, while 
because of the a coefficient in the second term in Eq. (41), 
ej. can be taken to the zeroth order. On the other hand, 
for less than 1 particle 

E = ae?. (45) 

Therefore, for the derivative jump at IV = 1 we have 

A E’ = - 4 (46) 


B. Singlet state with the number of particles between 1 

and 2 C. Triplet state with between 1 and 2 particles 


Let the spin-up state 0-|-(r) be fully occupied while 
the spin-down state </q(r, t) has the occupation a. Then 


Although it is also possible to analytically find n x -(r) 
in this case, the resulting expression is too lengthy and of 










6 


little use. On the contrary, the derivative discontinuity 
of energy at N = 1 can be easily evaluated. 

Let N = 1 + a. Then 


E = e 0 + ae 1 , (47) 

Pi+c*(r, r') = 0o(r)^ o (r') + a0i(r)0i(r'), (48) 

where <f>i and are the eigenfunction and the eigenen- 
ergy of the two lowest states in the potential v°^(r). 
Assuming a to be small, we can write by the perturba¬ 
tion theory to the first order in a 

E = e° 0 +ae° 1 + j Av(r)$ 2 (r)dr, (49) 


where <^(r) and e? denote the eigenfunctions and 
eigenenergies at a = 0, and Av(r) is the change of u(r) 
to the first order in a. From Eq. (16) we have 


In the case of the singlet state with between 1 and 2 
particles, (54) yields 


E = ef + aej. - J (f%(r)vf(r)dr 
-a j <f>l(r)vi(r)dr + a j 
Using the fact that 


^(r)^(ri) 

I r r 11 


(55) 


dr dri. 


5e c r 

Sv(r, a) 


</£( r ) 


and equating to zero the functional derivatives of Eq. (55) 
with respect to v(r,a), we have 


¥ 2 t (r') 

^(r,t) 



dr' = 0, 


(56) 


J Av(r)(j)Q 2 (r)dr + a J v° (r)^ 2 (r)dr+ 
a f <fto( r )^i( r )<fto( r i)0i( r i) ~ <^o 2 ( r )^i 2 ( ri ) drdr ^ = q 


T-ri 


(50) 


where u°(r) is v(r) at a = 0. The latter, however, is 
zero due to results in Sec. Ill A and the continuity of the 
potential at integer N for a fully spin-polarized system 
(see the end of Sec. II). Therefore, Eq. (49) yields 

E [+o = 4~ 

^o( r )^( r K( r i)^i(ri)-^ 2 (r)^ 2 (ri), , ( 51 ) 


drdri, 


J l r r iI 

and with the account of the results in Sec. Ill A 


r ') 


5( r '4) - 


0?(r") 


J 

from which we conclude that 


-dr" 


[ ^?(r') 

«(r.t) = a J y—^ dr ' + c ^ 

f <bl(r') 

5 M) = J y—^\ dv ' + c A 


dr' = 0, (57) 

(58) 

(59) 


with arbitrary c a . It is convenient to set these constants 
to zero, which makes the potentials zero at infinity and 
gives for the energy 


r)^(r') , 

E = e + + ae± — a —;-- —drdr . 

T * J |r-r'| 


(60) 


A E' = e\- 




t>°o( r)<A? 


( r i) 


r - r i 


drdri. 


(52) 


D. Comparison with Optimized Effective Potential 

In the case of OEP, we minimize the energy of the 
mixture 

E = (l-a)($iv|77Ar| < i)Ar)+Q(4>jv+i|7fAr+i|4 , Ar-i-i). (53) 

Subtracting and adding the effective Hamiltonian, by the 
same way as in Sec. II we arrive at 

N 

E = ^2,ti + ae N+1 - / riN +a (x)n(x)dx+ 

2=1 

1 r n, +a (xW(x) ,_1 /•|p jv + a (x,xQ| 2 

2 J I r r' | 2 J |r-r'| 

(54) 


Comparing Eqs. (39), (40), (36), and (41) with (58), 
(59), and (60) we observe two differences between LHF 
and OEP theory for a singlet with 1 < N < 2: (i) Differ¬ 
ent occupations of HOMO used in the calculation of the 
effective potential: It is the physical occupation a with 
OEP, while it is the re-normalized one, /3 of Eq. (37), 
with LHF; (ii) With OEP, the potentials are defined up 
to two arbitrary constants even at a fractional number of 
particles. With LHF, for the fractional case, the poten¬ 
tials are uniquely defined. Furthermore, since at integer 
N (3 = a, ior N = 1 and N = 2, LHF and OEP or¬ 
bitals and the energy coincide, while they are, generally 
speaking, different for 1 < N < 2. 25 We note, that with 
the method of Ref. 22, the energy with OEP can be also 
made the sum of the orbital energies, while the constants 
in the potential become fixed. 

In Fig. 2, we plot the energy of He ion versus the num¬ 
ber of electrons with the use of the LHF and OEP poten¬ 
tials, the differences between the two calculations being 
not discernible in the scale of the plot. Since at N = 1 
this is a fully spin-polarized case, the exchange potential, 
presented in the right panel, does not have a jump when 
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FIG. 2. Left: Total energy of an ion with the helium nu¬ 
cleus versus the number of electrons obtained with the LHF 
potential, and with OEP, the latter equivalent to HF for a 
singlet with 1 < N < 2. At N < 1, LHF, OEP, and HF coin¬ 
cide and are exact, the energy in that range plotted with the 
black straight line. Right: Spin-up exchange LHF potential 
at some numbers of particles close to 1. The potential does 
not experience a jump when N increases through 1 (see text). 


FIG. 3. Total energy of an ion with the beryllium nucleus 
versus the number of electrons, obtained with the potential 
of Eqs. (15) and (21), with the LDA xc functional, and with 
HF. 


close to HF, so in this respect LHF potential produces 
results very close to both HF and OEP not only for 
integer 18-20 , but for fractional number of particles too. 


the particles number changes through the integer value 
1, which is in accordance with the results of Sec. II. 

It is conceptually important that the scheme of the 
optimized time propagation, with the independent vari¬ 
ation of the orbitals instead of the potential, reproduces 
exactly the HF equations 26 . However, as we have seen, 
this scheme with the variation of the potential only, does 
not lead to OEP, but it rather leads to LHF. 


IV. NUMERICAL RESULTS AND DISCUSSION 

We have carried out self-consistent calculations of the 
total energy of beryllium and magnesium ions as a func¬ 
tion of electrons number with the use of the potential 
given by Eqs. (15) and (21). The specific ranges of the 
particle number variation have been chosen to keep the 
spherical symmetry. In both cases we do not extend the 
number of electrons above that of a neutral atom, since 
neither HF nor LHF support Be - and Mg - ions. 

Our results are presented in Figs. 3 and 4. For com¬ 
parison, we use the HF calculation performed with the 
use of Psi4 27 and NWChem 28 with aug-cc-pVQZ basis 
set, the latter two packages producing practically iden¬ 
tical results. For Be, we also plot the local-density ap¬ 
proximation (LDA) values 29 obtained with the correla¬ 
tion functional of Ref. 30. 

Our results, shown in Figs. 3 and 4, are essentially 
identical to the HF ones and are obtained with much less 
computational effort. On the other hand, it is known 14 
that in the present context, OEP also gives results very 
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FIG. 4. The same as Fig. 3, but for Mg ion. 

To explicitly demonstrate the jump in the effective po¬ 
tential, in Fig. 5 we plot the LHF potentials v x (r,a) for 
the Be ion with the number of particles 2.8 (lVj-=1.8, 
N±= 1), 2.9 (IV t =1.9, N±= 1), 3.1 (JV t =2, JV*=1.1), and 
3.2 (IVj-=2, IVj.=1.2). The potential varies smoothly and 
slowly from N = 2.8 to IV = 2.9 and from N = 3.1 
to N = 3.2, but it experiences a jump when the parti¬ 
cle number N passes through the integer value 3. The 
jumps of the potentials are plotted with the black dashed 
lines. It can be seen that they are constant in r in the re¬ 
gions, where the corresponding particle-densities are non¬ 
vanishing. It can be also verified that in these regions the 
jumps satisfy the relation 2Av x ^ + = 0, as required 

















r (a.u.) r (a.u.) 

FIG. 5. Spin-up (left) and spin-down (right) potential, calcu¬ 
lated by Eqs. (15) and (21), for Be ion for electron numbers 
2.8, 2.9, 3.1, and 3.2 as a function of the distance from the 
nucleus. The jumps of the potentials Av x (r, 4) and - A v x (r, 4.) 
between the number of particles 2.9 and 3.1 are plotted with 
black dashed lines. The solid black lines are particle-densities 
of spin-up and spin-down electrons at IV = 3. 


by Eqs. (30). In the asymptotic region, the jump for the 
spin-down exchange potential cannot be constant in r 
since, according to Eqs. (24), Vx 0+ (r, j,) ~ — 1/r+consti 
and Vx + ° + (r, 4) ~ const 2 . Further particulars of these 
calculations are presented in the Appendix C. 


of the results from the Hartree-Fock theory. Unlike the 
optimized effective potential, the LHF potential relies on 
the occupied states only and is, therefore, comparatively 
efficient and easy to implement. Surprisingly, we have 
found that a rigorous mathematical formulation of the 
LHF potential for fractional HOMO occupation dictates 
that the effective potential must be calculated with a 
renormalized value of the fractional occupation, /3 ^ a, 
where a is the physical occupation. Naturally, no renor¬ 
malization is needed at integer particle numbers. As a 
byproduct, our analysis shows the dangers lurking in an 
uncritical fractional filling of the HOMO, otherwise us¬ 
ing a theory derived for systems with integral number of 
particles. We have shown that with LHF potential the 
theory can be advanced much further in the explicit an¬ 
alytical form than it is possible with other methods. As 
a result, a deeper insight in the general aspects of the 
fractional occupation numbers theory becomes possible. 
Finally, we expect that our findings will boost the use 
of LHF potential in a broad area of applications where a 
consistent treatment of the fundamental gap is necessary. 
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V. CONCLUSIONS 


Appendix A: Details of the derivation of equations in Sec. II 

Introducing the reduced density-matrices 


We have demonstrated that the effective potential of 
Della Sala and Gorling, extended to the fractional num¬ 
ber of partcles, supports the derivative discontinuity, the Nr \ AH f ,, 2 , , 

Pk (xi,-X fe ) = r N _ k y J |$iv( x l--- x iv)| dXfc+i-dXjv, 


latter being a requirement of the exact many-body the¬ 
ory and the exact DFT. The quality of the results ob¬ 
tained from this effective potential is the same as that we can evaluate 

I 


($n\pn{*-)Un\®n) = 


P2( x > x ') . / , 1 f P3 ( X ’ X '> X ") , / 


. -dx + 
r-r' 2 


I r 


dxdx , 


($iv|/3jv( x )Vjv|$iv) = v(x)n N (x) + J p% (x, x')v(x')dx'. 
Therefore, Eq. (14) can be written as 

(1 — ^)Fjv(x) + /3-Fzv+i ( x ) = 0, 

where 


Fn(x) = J 


. - dx 1 , 
r-r' 2 


P2 ^ X ’ X ^ 1 1 J P3 ||, X ^ r/ X ^ dx'dx" - v(x)n N (x) ~ J P2 ( x > *')v(x)dx'. 


(Al) 

(A2) 

(A3) 

(A4) 

(A5) 


r 


To make further progress, we notice that the functions where 'L(x) are the standard field operators. Since the 
p of Eq. (Al) can be expressed as state |$at) is described by a single Slater determinant, 

(xi, ...x fc ) = (4> iv |4 ,t (xi)...^ t (x fc )4'(x fc )...4'(xi)|$ A r), 

(A6) 
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a straightforward application of Wick’s theorem 31 en- where 
ables us to express p £ as a product of densities and one- 
particle density matrices. In particular, for p^ and p% 
we obtain 32 


N 


/( x > x ') = J 2 <l>i( x )<f>i( x ') 


(A9) 


»=1 


P 2 (x, x') = n N (x)n N (x') - |p JV (x,x')| 2 , (A7) idempotent, i.e., 


is the one-particle density matrix. We note that the one- 
particle density matrix of a Slater determinant state is 


I p N (x, x")p N (x", x')dx" = p N (x, x') . (A10) 


Let us introduce the notations 

iV+7 = n 


(x,x',x") = p N {x,x')p N (x',x")p N (x",x)+ 
p N (x, x")p jv (x", x')p N (x', x) + n N (x)n N (x')n N (x' 1 )-- 

n JV (xO|p JV (x,x")| 2 -^(x")|p JV (xy)| 2 -n^(x)|p^(x',x") | | 2 , 

-—- 


p^ = (l-7)^+7^, 


(All) 


where 0 < 7 < 1. Then it follows from Eqs. (A4) and 
(A5) that 


P2 +P ( X ^') ± J ! 1 f P3 +P ( x i x ', x ") 


N+t3, 


dx! dx” — n(x)n JV+, 3 (x) — J p^ + ^ (x,x')v(x')dx' = 0 . (A 12 ) 


The key point in our derivation is that Eqs. (A7)-(A8) 
remain valid with p^ +1 in place of p^ on the left-hand 
side and p N+1 and n N+1 in place of p N and n N , respec¬ 
tively, on the right-hand side. The validity of the latter 
statement follows from the generalized Wick’s theorem 
proven in the Appendix B, or can be alternatively veri¬ 
fied by the direct substitution. This leads us immediately 
to Eqs. (15). Furthermore, since 

jF]\r(x)dxJ ^Y^——j^-dxdx' — Njv(x)n N (x)dx, 

(A13) 

integrating Eq. (A4), we have 

N{1-P)G n + {N + 1)PG n+1 = 0, (A14) 


or 

(1 — cx)Gn T ttGjv+i = 0 , (A15) 

where 

Gn = Jv(x)n N (x)dx ^y——■^-ypdxdx'. (A16) 

Using again Eq. (A7), from Eq. (A15) we arrive at 
Eq. (16). 

Finally, writing Eq. (15) with the explicit notations 
for the space and spin coordinates and noting that the 
density-matrix is diagonal in spin coordinates, we have 


[vx +P (r, <r) + G N+ p\ n N+f} (r , a) = J v^ +f} (r x , cr) - 

f Pa +/3 ( r , *i)Pa +/ 3 (ri, r 2 )p% +/ 3 (r 2 , r) 


r-n 


\Pa +0 ( r i r i)\ 2 dri 


(A17) 


In — r 2 


dridr2. 


If our system is not fully spin-polarized, then, for at least 
one spin direction (let us denote it with <5), the number of 
particles is integral and non-zero. In this case Eq. (A17) 
can be simplified. Integrating Eq. (A17) over the space 
coordinate r at a = S and taking account of the idem- 
potency of the density-matrix for the integral number of 
particles (p 2 = ps), we conclude that Gjv+p = 0 , which 
finishes the proof of the properties of the solution (15)- 


(16). 

Appendix B: Generalized Wick’s theorem 

Consider the expectation value of a product of cre¬ 
ation and destruction operators ABC....XYZ in a single 
Slater determinant state of N particles, denoted by |1V). 
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According to the standard Wick’s theorem 31 
(N\ABC....XYZ\N) = (N\AB\N)(N\CD\N)...(N\XZ\N) 
+ all possible pairing schemes, (Bl) 

where each pairing scene carries a sign plus or minus 
according to the parity of the number of interchanges of 
fermion operators that are needed to go from the original 
arrangement of operators to the paired one. There is 
no loss of generality in assuming that the creation and 
destruction operators ABC... refer to the same set of 
single particle states out of which we have selected the 
states cf>i, ..<f>jy that are occupied in | N) 

Consider now the state \N + 1} which is also a sin¬ 
gle determinantal state and differs from \N) only by the 
addition of one particle to single-particle state (f>N+i or¬ 
thogonal to 4>i,..(f>N- Again, according to the standard 
Wick’s theorem, the expectation value of ABC....XYZ 
in this state is 

(N+1\AbC...XYZ\N+1) = (N+1\Ab\N+1)... 

...(7V-|-l|AZ|A r -(-l}-l-all possible pairing schemes 

(B2) 

Combining these two expressions we see that the expecta¬ 
tion value of ABC....XYZ in the ensemble 7 |At + l}(A r + 

1| + (1_ 7 )|7V)(7V| is 

(1 - j)(N\Ab\N)...(N\XZ\N) 

+ 7 (AT + 1\Ab\N + 1 )...(N + 1\XZ\N + 1} 

+ all possible pairing schemes (B3) 

The generalized Wick’s theorem states that the above 
expression is equivalent to 

[(1 - ^)(N\Ab\N) + 7 (TV + 1\Ab\N + 1 )] 

...[(1 - i)(N\XZ\N) + 7 (N + l\XZ\N + 1}] 

+ all possible pairing schemes (B4) 

i.e., the standard sum of products of averages calculated, 
however, in the fractional ensemble. To prove the point 
we rewrite the last expression as 

[(A7| AB|AT) + 7 ((AT + 1|AB|AT + 1} - (A/"|AB|AT))] 
,..[(N\XZ\N) + 7 ((AT + 1\XZ\N + 1)(N\XZ\N))} 

+ all possible pairing schemes (B5) 

We need to show that, in spite of appearance, the above 
expression is actually linear in 7 . If this is the case, 
then it must necessarily coincide with the linear expres¬ 
sion (B3), since it obviously agrees with it when 7 = 0 
or 7 = 1. To prove that (B5) is linear in 7 we note that 
the quantity 7 ((AT + 1\Ab\N + 1 ) — (ARABIA/’)) for any 
pair operators A and B is either 0 or 7 . The latter case 
occurs when A and B happen to be, respectively, the cre¬ 
ation and the destruction operator of the state 4>n+i- It 
is also clear that such a pair of creation and destruction 
operators can appear at most once in the expression of 
ABC....XYZ. Therefore in each pairing scheme, there 
is at most one term proportional to 7 , and the whole 
expression is linear in 7 as we wanted to prove. 


Appendix C: Further particulars of the calculations 

In Fig. 6 we show two lowest spin-orbital energies of 
the Be ion versus the number of electrons. We also plot 
constants in the asymptotic behaviour of the exchange 
potential of Eqs. (24). No shift is given to the orbital 
energies, so the total energy is the sum of the latter 
[Eq. (12)]). While there are discontinuities in the orbital 
energies at integer number of particles, they compensate 
each other in the sum, the total energy being continuous, 
although having the derivative discontinuity (see Fig. 3). 
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FIG. 6. Orbital energies of an ion with the beryllium nucleus 
versus the number of electrons obtained with the potential 
of Eqs. (15) and (21). Open, semi-open, and solid symbols 
are energies of empty, partially occupied, and occupied or¬ 
bitals, respectively. The step lines show the constants in the 
asymptotic behaviour of the potential [Eqs. (24)]. 
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